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Abstract: In this paper, we solve the problem of recovering the motion of an observer 
relative to a surface which can be locally approximated by a quadratic patch directly 
from image brightness values. We do not compute the optical flow as an intermediate 
step. We use the coefficients of the Taylor series expansion of the intensity function in 
two frames to determine 15 intermediate parameters, termed the essential parameters, 
from a set of linear equations. We then solve analytically for the motion and structure 
parameters from a set of nonlinear equations in terms of these intermediate parameters. 
We show that the solution is always unique, unlike some earlier results that reported 
two-fold ambiguities in some special cases. 
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The problem of recoverying rigid body motion and surface structure uniquely from image 
data has been the topic of many research papers in the areaof machine vision [1-25] . Many 
approaches based on matching feature points [5,11,21,22] tracking contours [9j, and using 
velocity flow field [1,3,4,10,12,18,19,23-25], texture [2], or image intensity gradients [14- 
17] have been proposed in the literature, but not many have addressed the important 
issue of uniqueness. 

Feature point matching schemes require the detection of local brightness patterns 
that are likely to be found in consecutive images. A correspondence problem between 
the elements in successive 2-D images is first established and then the 3-D motion and 
spatial configuration of these isolated features is recovered. The minimum number of 
points required to recover the 3-D motion uniquely depends on the number of image 
frames. With 2 frames, in most cases, a minimum of 5 points results in a unique solution 
from a set of nonlinear equations. Using 8 points, however, as in algorithms proposed 
by Longuet-Higgins [11], Tsai and Huang [22], Buxton et al. [5], one only solves linear 
equations. In any case, these methods require feature detection and matching that fail 
to give reliable results when the object in view is smooth with no well-defined features. 
Further, since these methods use information only from a small portion of the image, 
they are noise sensitive. When nearby feature points are selected, these methods become 
even more sensitive to small amounts of error in the data. 

For smooth curved surfaces, Longuet-Higgins and Prazdny [10] suggested a method 
that uses the optical flow and its first and second derivatives at a single point. They 
reduced the problem to that of solving a cubic equation and concluded that, in general, 
three solutions are feasible. Later, Waxman and Ullman [23] developed the method 
proposed in [11] into an algorithm for recovering the structure and motion parameters 
from a set of nonlinear equations. They had to treat many special cases, and uniqueness 
results were shown through numerical examples. Waxman et al. [25] found a closed form 
solution to the original formulation in [23]. Apparently however, in their formulation, they 
are not able to resolve the two-fold ambiguities that arise in two special cases. When 
the observer translates parallel to the line of sight they cannot recover the curvature 
parameters. The result is a two-fold ambiguity similar to the one associated with planar 
surfaces. Also when there is no (component of) motion in the direction of zero surface 
gradient (labeled as structure motion coincidence in Waxman and Ullman [23]), they 
cannot distinguish between the two possible solutions. In any case, these methods are 
very noise sensitive since second order derivatives of errorful optical flow data are used. 
More robust algorithms that use the information from the whole region of the image 
plane have been suggested [1,3,4], but one still has to compute a dense flow field from a 
sequence of images. 

The flow field based approaches assume that a reasonable estimate of the optical flow 
is available. In general, the computation of the local velocity field exploits a constraint 
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equation between the local intensity changes and the two components of the optical flow. 
This only gives the component of the flow in the direction of the intensity gradient. 
To compute the full flow field, one needs additional constraints such as the heuristic 
assumption that the flow field is locally smooth [7,8!. In many cases, this leads to optical 
flow fields that are not consistent with the true velocity field. Since, in general, velocity 
based approaches are not robust, i.e., the solution may change drastically with only a 
small amount of noise in the data, one can argue that approaches that use an optical flow 
field that is computed through heuristic assumptions are apt to fail. 

In this paper, we present a method that uses the intensity values in a sequence of 
images directly, in order to recover the motion parameters as well as the local structure 
of the surface patches on the object. Since we do not compute the optical flow, we do 
not need to make any heuristic assumptions. We assume that the surface is textured 
or has surface markings, and that it is smooth so that we can approximate it in a local 
region by a quadratic patch. We give a closed form solution for the motion and surface 
parameters, and show that the solution is always unique. 

In a paper in press, we have extended our methodology presented in this paper to 
recover the motion of any smooth surface that can be presented by an n^^ order Taylor 
series expansion (with no restriction on n) uniquely. We show uniqueness and give closed 
form solution for the motion and surface parameters. 

2. Preliminaries 

We first recall some details about perspective projection, rigid body motion, the motion 
field and optical flow, and the brightness change constraint equation. This we do using 
the tensor notation (summation over repeated indices) and the following convention. A 
2-D vector v will be referred to by v,, and a 3-D vector u by {wi,W3}, where Wi is the 
component of u) parallel to the image plane. 

2.1. Perspective Projection 

Let the center of projection be at the origin of a Cartesian coordinate system. Without 
loss of generality we assume that the effective focal length is unity. The image is formed 
on the plane 2 = 1. Therefore, the optical axis lies along the 2-axis. Let {Xi, Z} be a 
point in the scene. Its projection in the image is {it}, where: 

' z 



2.2. Rigid Body Motion 

In the case of the observer moving relative to a rigid environment with translational 
velocity {<i,t3} and rotational velocity {«;,, wa}, we find that the motion of a point in 
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the environment relative to the observer is given by: 

—r- = -€ij{ZWj - AjW-i] - tt, '17 "^ ^ij^iWj - 13, 

where Cjy, the anti-symmetric permutation tensor, is given by: 

1 



''^ \ -1 



2.3. Image Motion Field and Optical Flow 

The projection of the 3-D velocity field, induced by a rigid body motion, onto the image 
plane is known as the image motion field. It is a purely geometric concept that is uniquely 
defined when the surface and the observer motion are given. It can be determined by 
differentiating the perspective projection relationship (X, = x^Z) with respect to time: 

dx{ 1 .dXi dZ 

Substituting for dXijdt and dZ/dt, and simplifying the results, we obtain: 

dx 1 

-jr = -^kji^ik + XiXk)wj + Xj€ijW3 + -^{xih - ti), 

where ^,y is the Kronecker delta symbol (it equals 1 when i = j and zero otherwise). It is 
important to remember that there is an inherent ambiguity here, since the same motion 
field results when distance and the translational velocity are multiplied by an arbitrary 
constant. This can be seen easily from the above equation since the same image plane 
velocity is obtained if one multiplies both Z and {ti,tz} by some constant. So we conclude 
that the surface structure and the translation motion parameters can be recovered only 
up to a scale factor. 

The optical flow is the apparent motion of the brightness patterns. An optical flow 
field is a vector field that shows how the brightness patterns at one instance of time can 
be transformed into that at the next instance of time. It is something we hope to estimate 
from an image sequence, but it is not uniquely defined. Additional assumptions must be 
made to guarantee a unique result. 

The distinction between the motion field and the optical flow is an important one 
that is rarely made. To recover motion, we need the motion field, but the best we can 
do is to estimate an optical fiow field from an image sequence. In extreme cases, the two 
will be very different: consider, for example, a perfectly smooth and uniform rotating 
billiard ball, or the motion of the shadows cast on a stationary object due to a moving 
light source. Under favorable circumstances the optical flow is identical to the motion 
field. Here, we assume that the motion flow field and the optical flow are the same. 
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2.4. Brightness Change Constraint Equation 

The brightness of the image of a particular patch of a surface depends on many factors. 

It may for example vary with the orientation of the patch. In many cases, however, it 

remains at least approximately constant as the surface moves in the environment. This 

is a good assumption when the surface has strong texture or surface markings, and is 

commonly used in work on optical flow. If we assume that the image brightness of a 

patch is indeed constant, we have 

dE 

or, using the chain rule: 

dE_dx^ a£ _ 

dxi dt dt 

where Ex, = dEjdx^, and Et = dEjdt are the spacial and temporal derivatives of the 
the image brightness function. Substituting for dxi/dt we get: 

Et - {Sik + XiXk)tkjEx,vjj + Xje,jEx,W3 + -{xiExM - ExM) = 0- 
If we let d = 1/Z, and define: 

Si = —Exi, Sz = XiExi, 
V, = (.ij{xjSz - Sj), Vz - -tijSiXj, 

then: 

Et + Vm + Vz^z + d[S,U + 53*3) = 0. 

This is the brightness change constraint equation in the case of rigid body motion. 

3. Surface Representation 

In a small neighborhood, a smooth surface Z{Xi) can be represented by its Taylor series 
expansion about a reference point. Without loss of generality, let us consider a neigh- 
borhood of the line of sight, denoted by Tz, and choose the point of fixation, O, located 
at a distance Z„ from the viewer as the reference point. Clearly, the origin of the image 
plane is the image of the point of fixation. Then we can write: 

Z(X.) = 2„ + (^)„X..l(^)„X,X,+.... 

If we neglect the higher order terms, we are left with a quadratic approximation to 
the surface (see Figure la). While this representation is globally valid for quadratic 
surfaces, it may only be locally valid for other smooth surfaces. We can express the series 
expansion in terms of the image coordinates, x,, by employing the perspective projection 
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Figure 1: (a) A Quadratic Approximation to the Surface Patdi at P 
(b) Image of a Surface Patch and Regions of Amdytkity 



relationship (X, = XjZ). We consider the series expansion of d = 1/Z instead, since 
the brightness change equation is in terms of d. The surface is then represented by: 

where it can be shown (see appendix 3) that: 



' dxi 



z,:dXi 



' dxidxj 



' dXidXj 



Since the surface parameters and the translation motion can only be recovered up to a 
scale factor, we can arbitrarily fix one of the coefficients in the above series expansion. 
We choose to set d„ = l/Z,, = 1. Therefore: 

rf(z,) = 1 + dji, + -dijXiXj + ...., 

it 
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where 

In this paper, we will assume that the above series expansion, up to second order terms, 
is a good approximation to 1/Z in Tz- 

4. Essential Parameters for a Quadratic Patch 

Let Az denote the image of the surface patch Tz. We assume that the patch has texture 
or surface markings. More precisely, we assume that in small regions in Az corresponding 
to areas of texture, the intensity variation is smooth so that locally, we can approximate 
the intensity function by a Taylor series expansion. Equivalently, the image intensity 
function is assumed analytic in small areas that we refer to as regions of analyticity. A 
region of analyticity is contained within its boundary of analyticity (a similar term has 
been used by Waxman and Ullman [23] in a different context). This is illustrated in 
Figure la, where each bounded area in the image represents a region of analyticity. 

Now, consider the Taylor series expansion of the image intensity function, E, and its 
temporal derivative, Et, about a point p in an analytic region Api 

E = E„ + EiSxi + -EijSxiSxj + -EijkSxtSxjSxk + —EijkiSxiSxjSxkSxi + ..., 

Et =T„+ Tibxi + -TijSxiSxj + -TijkSxiSxjSxk + —TijkiSxiSxj6xk6xi + ..., 

where 

6xi = X, - xf . 

If we define: 

E„ = E„ - E,x\ + \e,,x^ - \E,,kxlx]xl + ^^..xfxjx^.xf - ..., 

Ei = E, - E,^x) + \E,,kx]xl - \E,,kix]xlx\ + ..., 
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Eij - Eij - Eijkxl + -E^jkixlx^i - ..., E^jk = Eijk - Ei^kix'^^ + ..., Ei^ki = E^jki 



J,, - J,, Jt^t -^ 2 '-^ ' > Q-^^J^x^XjXj^ + ^^i^jklx^x-Xf^Xi ..., 

Ji — J^i ~ -^ijXj + l^^ijkXjX^. — -TijklXjXf,!^ + ..., 

Tij = Tij - TijkXf^ + -T^jkiXkXi - ..., Tijk = Tijk - Tijkix^i + ..., T^jki = ^t;fc/ 
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then for the neighborhood Ap, we can express the Taylor series expansion of the intensity 
function and its temporal derivative, about the origin as: 

E = E„+ E,x^ + -E^jXiXj + -E^jkXiXjXk + —E.jki^iXjXkXi + ..., 

Et = T„-^ TiXi + -T.jXiXj + -TijkXiXjXk + —Tijkl^i^j3:kXi + .... 

Each region of analyticity will give rise to an independent Taylor series expansion for 
the intensity function and its temporal derivative. We will show that if the patch is 
strongly textured such that higher order coefficients of the intensity function, as in the 
above series expansions, exist and can be computed realiably, then we can compute the 
appropriate parameters from the information in only one region of analyticity. For rea- 
sons that become apparent, if we restrict ourselves to second order terms of the Taylor 
series expansion, then we need at least two such regions to recover the surface and mo- 
tion parameters. In most cases, the texture varies more rapidly than the surface, and 
therefore, several such regions exist within each quadratic surface patch. This is a rea- 
sonable assumption for most textured surfaces. Similar assumptions are generally used 
to determine structure from motion. When several regions exist, then a least-squares 
formulation can be implemented to exploit the information from every region. 
The spatial image intensity gradient is given by: 

Exi = Ei + EijXj -(- -EijkXjXk + -EijkiXjXkXi + —EijkimXjXkXiXm + ■•■ • 

Using the expression for the image intensity derivatives in the definition of {Si,Sz} and 
{^^^3} we arrive at: 

Si = -{Ei + EjjXj + -EijkXjXk + -E^jklXjXkXi + —E^jklmXjXkXiXm + ....), 

S3 = EiXi + EijXiXj + -EijkXiXjXk + -EijkiXiXjXkXi + ...., 
Vi = iijEj+eijEjkXk + UjEkXjXk + -(ijEjkiXkXi + tijEkiXjXkXi+ 

-^ijEjklmXkXlXm + -eijEklmXjXkXiXm + —€ijEjklmnXkXlXmXn + ■•••, 
^3 - -{tijEjXi + eijEjkXiXk + -tijEjklXiXkXi + -tijEjklmXiXkXiXm + ....)• 

Substituting for {5,, 53}, {Vj, V3}, d, and Et in the brightness change constraint equation, 
and ignoring higher order terms, we get: 

1 1 1 

P„ + PiXx + -PtjXtXj + -p.jkXiXjXk + —pxjklXtXjXkXl = 0, 



g Direct Passive Navigation 

where 

Pi, = 7^, -t- nijEj, 

p, - T, + m^jE.^ + nijEj^, 

p,j = r,y + m^jkEk + 2m^kEk] + nikEktj, 

Pi]k = T^jk + '^ij-E'*: + ^^ijkEit + SmikEitj + miE^jk-, 

Ptjkl - Tijkl + ^rkijEkl + ^rrirmjEmkl + ^T^mtEmjkl + 'mmEmijkh 

and 

rrii - ejiWj - ti, 

rriij = Sijts - tijWz - ditj, rhtj = tsd^j, 

'^tjk = {^ij^l + djt3)Ski + {iki^k + dit3)6kj - tkd^j. 

This is the Taylor series expansion of the brightness change constraint equation for the 
motion of an observer with respect to a locally quadratic surface. We will term m^, 
m,j, fhij, and rriijk the essential parameters (these are similar to the ones we defined in 
an earlier paper [16] for the case of the motion of an observer with respect to a planar 
surface, in agreement with Tsai and Huang [21]) since as we see later, these are essential 
to recovering the motion and surface parameters uniquely. We will first compute these 
parameters using the image intensity gradients in local neighborhoods, and then use them 
to recover the unknowns {{ti,tz} ,{wi,wz}, and {di,dij}). 

We should remind the reader that if we do not restrict ourselves to quadratic patches, 
then the expressions for pijk, PijkU etc, will include the higher order terms of the surface 
[dijk-, dijki, etc). In a paper in press, we have extended the results presented here to 
derive a unique solution in closed form for the general case of smooth surfaces that are 
represented by an n*^ order Taylor series expansion (with no restriction on n) . 

4.1. Constraint Set 

For the brightness change constraint equation to be satisfied at every point in Az, all of 
the coefficients of its Taylor series expansion must vanish. That is, 

r„ + mjEj = 0, 

Tj + rUijEj + rUjEji = 0, 

Tij + rriijkEk + {m.kEkj + rn^kEki] + rrikEkij = 0, 

T^jk + rh.jEk + {mijkEu + mi^Eij + mi.jEik) + {mikEuj + mijEuk + mi,Eijk)+rniEujk = 0, 
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Tijkl + [rhijEki -r fh^kEji + rh^iE^k ^ rhji,E,i ^ rhjiE^k ~ rkkiE^j)- 

[mmijEmkl + mrmkEmjl ^ mmilEmjk + rrimjkEmU ^ ^mjlEmik ^ ^mklEmi]^ 
{rrimiEmjkl + '^mjEmikl + ^mkErmjl + ^mlEmijk) + rnmEmijkl ~ 0- 

These are 15 linear constraint equations (note that the last 3 tensor equation only give us 
3, 4, and 5 constraints because of the symmetry in the indeces) in terms of the 15 essential 
parameters (two mj's, three m^y's, four mj/s, and six m^jfe's). We will refer to them as a 
constraint set. Each region of analyticity will give rise to a constraint set (as opposed to 
a brightness constraint equation for each point in the image). Our strategy is to, first, 
determine the 15 essential parameters from the image derivatives using a constraint set. 
If several such regions exist within the quadratic patch, then a least-squares formulation 
is used so that we can exploit the information from more areas of the image. We then 
compute the motion and surface parameters from the essential parameters. 

We remind the reader that there are only 11 surface and motion parameters to recover 
{{iiM}-, {«^t,«^3}) and {di,dij}), however, we are measuring 15 essential parameters. We 
will show how the extra information will allow us to resolve the resulting ambiguities that 
arise in some special cases. 

5. Recovering the Extended Essential Parameters 

Let p be a point in Az, and let: 

1 1 1 ^ 

Ex, = E^ + E.jXj + -E^jk^jXk + -E^jkl^jXkXi + —E,jklm^i^jXk^l^m + •••• 

Et = r„ + r,ii + -T^jX^Xj + -TyfcXjX^Xfc + —Tijki^^t^j^k^i + ••■• 

represent the Taylor series expansion of the intensity derivatives m Ap, a. region of ana- 
lyticity in Az around the point p. Theoretically, the coefficients of these polynomials, up 
to fifth order terms, can be determined using the intensity data in Ap through a least- 
square-errors fit if the patch is strongly textured. We will show, however, that this is not 
necessary for determining a unique solution for the essential parameters. 

Assuming these coefficients can be computed for at least one region of analyticity, Ap, 
we can solve for the 15 essential parameters, by elimination, from the 15 linear equations 
in the constraint set for region Ap-. 

T„ + m^Ej = 0, 

r, + ruijEj + rrijEji = 0, 

Trj + m.jkEk + [mrkEk] + rUj^Eki) + rUkEkij = 0, 

Tzjk + rh^jEk + {mijkEu + mi^kEij + mujEik) + [mikEuj + mijEuk + rnuEijk) + miE^jk = 0, 
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Tijki + irfixjEki + fh.kEji + rhaEjk + rh^kE^t - m^;^,* + fhuEij) 

+ {mmijEmkl + ^mikEmjl + T^mxlEmjk + rrimjkEmil + "^mj/^'mtA: + rUmklEmi]) 
+ {mmiEmjkl + rUmjEmtkl + mmkEmtjl + rnmlEmi]k) + mmEmijkl = 0- 

When several regions of analyticity exist on the surface patch, then a least-squares for- 
mulation should be implemented to exploit the information from every region. 

As we said above, it is not necessary to compute up to the fifth order terms of the 
intensity function. We can get an approximate solution from: 

T„ + rrijEj = 0, 

Ti + rriijEj + mjE^i = 0, 

Tij + rriijkEk + [rriikEk] + rrijkEki) + rnkEki, = 0, 

T,jk + fn,jEk + {mijkEu + mukEij + mujEik) + {mikEuj + mijEuk + mi.Eijk) = 0, 

{fh^jEki + fh^kEji + fh^iEjk + fhjkEti + mjiE^k + fhkiEij) + 

{rrimtjEmkl + rUmikEmjl + rrimtlEmjk + mmjkEmil + T^mjlEmtk + mmklEmij) = 0, 

that is, using only up to the third order terms of the intensity function. Similarly we 
may drop the third order terms and use the resulting equations to compute the essential 
parameters. In essence, we are compromising the accuracy of the solution of the essential 
parameters for two objectives. First, it is more practical to fit quadratic or cubic poly- 
nomials to intensity function in a small neighborhood, and secondly, less computation is 
required to do so. 

We will show that in all but two cases, m^, m,y, and m,jk are suflacient to derive 
uniquly a closed form solution for the motion and surface parameters. For these two cases, 
we have to resort to the constraints given by m,y to resolve the two-fold ambiguities in 
the solution. In other words, we can use the equations for m,, m,_,, and niijk to compute 
{wt,W3}, {t,,<3}, and {<i,,d,_,} uniquely except for two special cases, when these arise, 
then we use the equations for m^y to resolve the resulting two-fold ambiguities. 

Since the two special cases may occur rarely, we can consider a different strategy that 
does not require the computation of m,j. Since fhij only appears in the last two tensor 
equations, we consider only the following equations: 

T„ + rn^Ej = 0, 
Tj + rriijEj + nijEji = 0, 

Txj + m^jkEk + [m.kEk] + m^kEki) + rrikEktj = 0. 

We can use these equations to compute the remaining essential parameters, however, we 
only have 6 constraint equations, for each point p, in terms of the 12 unknown essential 
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parameters. So we need at least two points. Using two points, we can solve for m^ as 
follows: 



G.^^^ 



where 



Gm — 



'pi 



1 / El -E\ 



''' ^^kE]El\-El E\ 



and superscript p denotes the coefficients associated with point p. Using the solution for 
mi, we can now solve the second tensor equation for rriij. The solution is given by: 

m,, = -G,,{T^ + El^m,). 

Finally, using m, and rriij from above, the solution for rriijk can be expressed by: 

m^jk = -Gpk[^ + E^i^mu + Elrriji + miEl^). 

It is now clear why we need to consider two points from two different regions of analyticity. 
To be able to solve for the essential parameters, it is necessary and sufficient that: 

U,E}El^O. 

If the two points are chosen from the same region, the coefficients of the Taylor series 
expansion of the intensity function about each point will be the same, and therefore: 

U,ElE}^0. 

It is appropriate to point out what information is necessary to recover the essential 
parameters. It is necessary to compute the first order coefficients of the image function, 
however, the higher order terms only provide more accurate estimates of the unknown 
parameters. To show this, let us assume that only the first order terms are available. 
Then, we can still determine the essential parameters from: 

m, = -Gp,rf , m,, = -Gp;Tf , m,,k = 0. 

The terms we neglected could be viewed as the correction terms that come from the higher 
order information. We can still estimate the motion and surface parameters without 
them, even though these estimates may not be very accurate. Similarly, we can get an 
approximate solution for ruijk by neglecting the last term involving Euy In this case, we 
only have to compute up to the second order terms of the Taylor series expansions of the 
intensity function and its temporal derivative in each region of analyticity. In practice, 
this is desired not only from a computational consideration, but also with noisy images, 
one cannot determine these coefficients accurately. 

Image brightness values are corrupted with sensor noise and quantization. So it is not 
advisable to base a method on measurements at just a few points. Instead we propose a 
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least-square-errors formulation. We choose m,,, m,j, and m^jk '^hat minimize: 

{Tij + m,jkEk + {m^k6mj+rnjk6mi)Ekm + rniEuj) j dA, 

over the whole region Az. For an extremum of J, we must have: 

dJ dJ ^ dJ ^ 
= 0, = 0, and = 0. 



<9m, dm^j drriijk 

After differentiating the expression for J with respect to the essential parameters, and 
performing some algebraic manipulations, we arrive at: 

Ak"^k + ^ikl^kl + ^iklm^klm = -^t» 
Bijk'^k + Bijkirriki + Bijklm"^klm - -Bij, 
C,jkii^l + Cijkim^lm + Ckirriiji = -Cijk- 

where 

Ai = I {T„E, + T,Ea + 2Tr,E, + Ti,Euj)dA, 
Ja. 

Aik = f {E,Ek + E.jEkj + dk^Eij,qEi^q)dA, 

Am = / {EiE,k + 2EkjE,ij)dA, Aikim = I {EkEum)dA, 
Ja, Ja, 

B,, = / {T,E, + 2T,kEjk)dA, B^jk = f {Ek^Ej + 2E,iEkn)dA, 
J A, Ja, 

Bijki ~ I [kkEiEj + 26ikEimEjrn + '^EjkEii)dA, 
Ja, 

Bijkim = / ihiEji + 6iiEjk)EmdA, 
Ja. 

C^jk = f {Tr,Ek)dA, C,,kl = I {EkEuj)dA, 

J A, J A, 

Qjkim = / [StiEmjEk + SjiEmiEk)dA, Cki = / {EkEi]dA, 
Ja, Ja, 

These twelve linear equations can easily be solved for the optimum rrii, rriij, and m^jk 
by elimination. The derivation of the results is a tedious practice of elementary algebra 
and is omitted here (note that the last tensor equation only gives 6 independent scalar 
equations due to the symmetry in i and j). 
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So in summary, we have shown how to compute the essential parameters from just 
one point (using the fifth order terms of intensity function in the region of analyticity 
containing the point), or two points (using third order terms of intensity function in the 
region of analyticity containing the point). We also showed that approximate results may 
be obtained if one uses just the lower order terms. When more points (or correspondingly 
more regions of analyticity) are available, a least-squares formulation can be employed 
that exploits the information from larger portions of the image plane. 

In most cases, it will be enough to compute up to second order terms of the Taylor 
series expansion of the image function, however, as explained earlier, we have to resort 
to the information provided by the higher order terms to resolve the resulting two-fold 
ambiguities in the solution for two special cases. 

6. Recovering the Motion and Surface Parameters 

The essential parameters were defined in terms of some nonlinear functions of the motion 
and surface parameters, {^1,^3}) {w^t)W^3}f and {di,dij}: 

TUi = CjiWj - ti, 

rriij — dijtz — €ijW3 — ditj, fhij — t^dij. 
^ijk = {(ij^l + djtz)6ki + [tiiWi + ditz)Skj - tkdij. 

We will show that the equations for m^, rriij, rriijk can be used to recover the motion 
and surface parameters uniquely in all but two cases. When there is no component of 
motion parallel to the image plane (i,y (these are related to the curvature of the surface 
around the fixation point) cannot be recovered, and therefore, we have to contend with 
a two-fold ambiguity similar to the one associated with planar surfaces in motion, and 
when there is no motion in the direction of zero surface gradient, there is yet another 
two-fold ambiguity. It is in these cases that the equations for rhij come to rescue. As a 
result, we are able to recover the motion and surface parameters uniquely in every case. 
So we will not use the equations for fhij until the very end. We choose this strategy since 
to compute m^ , as we showed earlier, we need good estimates of the higher order terms 
in the Taylor series expansion of the image intensity function. In practice, these terms 
cannot be computed reliably from noisy images. 

It is possible to solve the above non-linear equations for the motion and surface 
parameters in several ways. Curiously, they all reduce to solving a cubic equation in 
terms of ti/t2 (or <2/*i) that was first derived by Longuet-Higgins and Prazdny [] in 
their work using optical flow (however, they never tried to solve it analytically, and so 
concluded that, in general, three solutions are possible). Once the cubic equation is 
solved, the unknown parameters can be determined easily. 

It was brought to our attention that Waxman et al. [25] had derived the same results 
from optical flow data using the formulation originally proposed in Waxman and Ullman 
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23;, and by exploiting 12 observables that are actually linear combinations of our 12 
essential parameters m,, m,j, and mi^jt), that are computed from the optical flow and its 
first and second derivatives. Unfortunately, their formulation, using the 12 observables, 
does not allow them to resolve the resulting two-fold ambiguities in the two special cases. 
Here, we present a procedure that uses the image coordinate transformation that 
Longuet-Higgins and Prazdny suggested. Before we proceed, we will redefine rriijk by 
subtracting twice the expression for rrii from the old definition of m^jk. The new equations 
will be given by: 

rrii = tjiWj — ti, 

rriij = bijtz — tijWz — ditj, rriij — Hdij. 

^ijk = {tj + djt3)6ki + {ti + ditz)6kj - tkd^j. 

An image coordinate transformation can be written as: 



, _ _ / cos — sin 6 \ 



We will explain how 6 is chosen later. Let m'-, m'--, and m^ ^^ denote the extended essential 
parameters in the new coordinate system. Then: 

m[ - Qiimi = il,ii{ejiWj - <;), 

m[ji^ = i^il^jm^kn^lmn = ^il^jm^kny{tm + dmt3)Kl + (</ + ditz)6nm - <nC^/m j • 

These can be written in the form: 

^'x]k = ^kx{t'j + 4*3) + ^kj[t[ + d%) - t[d\j, 
where, 

w'- = UilWi, t\ — Qiiti, rf, = Vlildi, rfjy = ^il^jm^lm- 

Since the rotation to the new coordinate system is about the z-axis, it is clear that tz = t'^, 
and a'3 — w'^. 

Let us write the equations for the transformed extended essential parameters in the 
component form: 

171*2 = u;J — t'2 
"^'21 = "^3 - 4^1 

"^222 ~ ^1^2 "*" °2^3) ~ ^2^22 



m\ = -w'2 ~ t\ 




m'l, = ^'3 - d\t\ 




m\2 — -w'^ — d\t'2 




m\,, = 2{t\+d\t'^)- 


- t'A, 
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mjj2 — "^2*^11 "^221 ~ ~^i"22 

"^121 ~ '2 "*" ^2^3 "~ ^1*^12 "^122 ~ ^1 "*" °1^3 ~" ^2^12 

The new coordinate system is chosen such that the Xj-axis (or alternatively, Xj-axis, 
due to symmetry in the above equations) is aligned with the component of the translation 
in the image plane so that t[ — (if we align the Xj-axis with <,, then <2 — 0). Note that 
6 is not known, since t^ is unknown. We will show how to compute analytically. 

Setting t\ = 0, we arrive at: 

m'j = —W2 m'2 = w[ — t'2 

mjj = <3 17122 ~ ^3 ~ "2^2 

'^la ^ -^'3 - ^'A "^21 = "^3 

"^111 ~ 2a^t^ 1^222 ~ ^'^2 "*" ^2^3) "" ^2^22 

mjj2 = ~^2^ii "^221 ~ ^ 



m 



121 



— tj ■+" "2^3 '^122 ~ "1^3 *2"'\2' 



We use the constraint 771221 = to solve for 6. Writing this in terms of the original 
parameters, we obtain: 

After substituting for Dt^'s, collecting terms and simplifying, we arrive at: 

— mii2tan^^ + (mm - 2mi22)tan^d — (m222 - 2mi2i) tan^ + m22i = 0, 

or in terms of the motion and surface parameters, 

^2^11 tan^d+ {2t2di2 - firfn) tan^ ^ + (<2«^22 -2<idi2)tan^ -<i<i22 =0, 

which can be factored into the following form: 

(<2 tan e - ti){dn tan^ 6 + 2di2 tan 9 + d22) == 0. 

The solutions of the above equation are given by: 



ti -<^12 ± Y <^12 ~ <^11^22 



tan = —, 

<2 Oil 

The cubic equation is trivially satisfied for two special cases: d^j = 0, that is, the surface 
is planar, or t, = 0. If the coefficients of the cubic equation do not vanish, we are 
guaranteed that: 

ti 7^ 0, and rf,j y^ 0. 

For now, we assume these conditions are satisfied. 

Note that only the first solution, tan^ = ti/t2, corresponds to the correct transfor- 
mation. Let us look at the other two solutions, i.e., the solutions of: 

dii tan^ + 2di2 tan ^ + ^22 = 0. 
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Remember that: 



and therefore, 



where 



d,j = -( 



d^Z 



dXidXj 



u 



{d\\ tan^^ + 2di2 tan^ + ^22) = -VtH^jVj, 



Hij - 



Vi — ( tan 6 1 ) , 
d^Z 



dXidXi 
d'^Z 



\ 



\dX,dX2 dxl ) 

As shown in appendix 1, if the surface has a positive gaussian curvature at O, that is O 
is an elliptical paraboloid point, then Hij is either positive or negative definite, and so: 

du tan^ 6 + 2di2 tan 6 + d22 - -ViHijVj ^ 0. 

Therefore, the solution for the transformation angle is unique . If the gaussian curvature 
at O is zero but dij / 0, or if O is a parabolic cylinder point, the above quadratic equation 
has two identical solutions given by: 

tan = — —— — — —— . 

"11 012 

Finally if the surface has a negative gaussian curvature at O [hyperbolic paraboloid point), 
then the quadratic equation has two distinct real solutions that are: 



tan^ = 



-dn ± \Jd\^ - dnd22 



d 



11 



So the solution for the transformation angle is not unique when point O is either a 
cylindrical paraboloid or a hyperbolic paraboloid point. We will show that, in most 
cases, it is possible to single out the correct solution, 6 = tan~^ ^1/^2, using the extra 
constraint equations. 

Once we solve for the transformation angle(s), we can determine the transformed 
parameters, m|, mjy, and m'-'f^, and determine the motion and surface parameters from: 



m, 



— w 



m 



m 



m 



m 



m 



m 



11 " 
I 

12 " 

111 
112 
121 



4- 



d'A 



— w 
-- 2d\t 
- -t' d' 

- Co T*" t*o C o 



2=<- t'2 



m 



22 



^3 ^22 



m'21 = w'^ 



"^222 ~ ^'^2 ^ ^^2^3) *2^22 
= 



m 



221 



'"122 ~ "k3 '^2'*12' 



1^2 


^-m\ 




t',- 


= "^121 




d\- 


Ttl^] + '^12 


"^121 




d'n 


m',91 





w'^ = m'21 




t'z = rn[, = 




4 = -^ 
"^121 




jt "^222 


- 2m', 2, 


*^22 - 


0'. 



6. Recovering the Motion and Surface Parameterf 

If tz = t'^ = m'ji =0, the solution is given by: 

w[ = m'2 + m'i2i 
t[ =0 

d'o = l 

d' = ^112 
11 m'i2i *' »T^i2i " '■"121 

If tz = t^ — m'lj # 0, we determine w'2, 1^3, and ^3 as before: 

w'2 = -m'l, ,11)3 = mji, ts = mil- 

Substituting for these into the remaining equations, we obtain: 

m22 = "^11 — "2*2' 
mi2 = — »^21 "~ ^1^2' 

mm = 2mii<ii, 

''^222 — 2(^2 + "^11*^2) ~ *2<'22' 
"^112 — ■"^2*^11' 
"^121 = *2 + "^11<^2' 
"^122 = "lijrfj — ^2^12. 

We can solve for dj from: 

A = ^. 

^ 2m'ii 
Now expressing dj in terms of t'^ we get: 

.; _ "^121 ~ ^2 

di - — ~i ' 

mil 

and substituting into the equation for m22, we arrive at: 

^1 — ^1 ^2("^121 ~ ^2) 
m22 - "^11 - , 

This simplifies into: 

t'i - m'uA + ^11(^11 - "^22) = 0- 
The two possible solution for t'^ are given by: 



*2 = 2 ("^121 ^ \/"^'ni - ^m'lJm'ii - m'22)), 
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or in terms of the motion and surface parameters: 

<2 = t'2, and f^o's- 

We now proceed to determine the solution for the remaining parameters in terms of the 
solutions given for ^'2: 

., "^112 j> 2m'i22 - m'lii , _ m'222 - 2m;2i 

4 -^4 2 

Now consider the extra equation we have not used yet, i.e., 

mi2 + Tn2i = ~ayt2. 

In most cases, of the two possible solutions of t'2 given earlier, only the corect one satisfies 
this constraint equation. So the correct solution for t'2 (that will be unique) is chosen 
and used to determine the other parameters. Similarly, in most cases, this equality will 
not hold for the wrong transformation angle 6. So the above constraint can be used to 
determine the correct solution for t'2 (from the dual set) and the proper transformation 
angle (if the cubic equation we derived earlier has multiple solutions). 

In the special case when d'^ = 0, we have: 

m'i2 + m'21 = 0. 

So the extra equation becomes independent of t'2, and cannot be used to determine the 
correct solution of t'2. However, the constraint equations: 

come to rescue here. Of the two solutions for ^2 that result in: 

., "^'112 ji _ 2m; 22 - ^111 ., _ >^222 - 2^21 

On = --77-"' °i2 - :^, ' "22- ^, 

4 ^h '2 

only one is consistant with the solution obtained from: 

,t _ ^11 .1 _ "^12 ./ _ "^22 

"11 T"' "12 — I ' '*22 — I • 

mil "^11 "^11 

We again emphasize that the constraints from mj_, are only used for removing the ambi- 
guities, and not to compute the unknown parameters for the reasons given earlier. 

If d' = 0, the ambiguity is resolved without the need for the constraints from mj^. 
This is because the dual solution for t'2 (4^3) will be zero, and we can discard it since we 
assumed that t'2 ^ (<', ^ and t'j ^ imply that t'2 ^ 0) . 
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Now we consider the case when the coefficients of the cubic equation vanish. Then, 
either f, = 0, or dij = 0. Again, we use the constraints: 

^!; = i'sd'i,, 

to distinguish between the two cases, and determine the solution. 

If m'-j = 0, then d'-^ = since we assumed that t'^^O (the case ^3 = was treated 
earlier). The remaining parameters can be determined using the method we proposed in 
[16], however, the two-fold ambiguity of planar surfaces is what we have to contend with. 

If mj^^ 7^ 0> then we know that d[^ ^ 0, and therefore, t[ = 0. In this case, the solution 
is given by: 

W2 — — m'j iWg = m'jj = —"^'12 

t'2 = f^ = m'u = m'22 



w[ = 1712 




t[ =0 




d'o = l 




^1 ^\ 1 


^11 



' _„' _-' ™' 



m,,, mm ., moo, m 



ftl - "Mil _ "Mil Jl _ ^"-222 



2mii 2m22 ^ 2m\i 2m\ 



221 



22 
J' _ "^12 _ ^12 J/ _ ^22 _ ^22 

"12 — / — / "'99 — / — * 

mil "^22 ^u '^22 

This is the only case we actually need m,j to estimate the surface parameters d^y 
So in summary, we proceed as follows: 

(1) We determine the solution(s) for the transformation angle from: 

-mii2 tan 6 + (mm - 2mi22) tan^ B - (2mi2i — ^222) tan 6 + m22i = 0. 

If all of the coefficients of the above equation vanish, then either d^j = rf| = 0, or 
ti = t[ = 0. We use the constraint equations: 

to distinguish between the two cases. If m' ■ = 0, we conclude that d[- = and 
therefore, the surface is planar. We can use our method for planar surfaces (see 
[16)) to compute the unknown parameters, we have to contend with the well- 
known two-fold ambiguity of planar surfaces. If m', 7^ 0, then the solution is 
given by: 

w[ = m2 W2 — —m[ itig = 171*21 ~ ~"^i2 

t\ =0 ^2=0 t'^ = m'u = m[ 



22 



d'= I d' = - "^^11 = -"^111 d' = -"^222 = _nh22. 

^ 1 2mii 2m22 ^ 2mii 2m22 

^/=^_^ fV -Till -T!ii ■> _ ih'22 _ m22 

11 ' — ' "12 — I — / "■99 — / — / 

mil ^22 "^11 "^22 ^n "^22 
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When the coefficients of the cubic do not vanish, we conclude that: 

U = t[ # 0, and d^j = d[^ ^ 0. 

For each solution of 0, we compute the extended essential parameters in the rotated 
coordinate frame. 

(2) We determine the solution for t'^ from t'^ = m'jj. If t'^ - 0, then the solution for 
the motion and surface parameters is unique and is given by: 



■w[ = m'j + mi2i 


w'2 = -m'j 


w'z = "^'21 


t\ =0 




tj = mj2i 


«3 = ^'11 = 




.^112 
"^121 


, m'21 + m\^ 
mi2i 

J/ "^122 
"^121 


"^121 

J/ "^222 ~ 2mi2i 
22 - ™' 

"^121 



If <3 7^ 0, we first find the two solutions for t'2 from: 

<2 = 2('"'^" ^ y^mi2i-4m'„(m'ii -m'jj)). 

If one of the solutions is zero, we discard it (resolvable ambiguous case), and 
proceed with the correct solution. If both solutions are nonzero, we pick the 
correct one that satisfies: 

mi2 + m2i = ~flit2- 
If this constraint is trivially satisfied, that is: 

mi2 + »T^21 - 0' 

we conclude that d\ — 0. Of the two solutions given by: 

W2 = -m'j w'^ = m'^i 

t\ = f'2, and 4<^ t'^ = m'u 

^ 2mii -^ m'u 

, _ 2m', 22 -^111 J/ _ _^222 -2m,2i 

"12 ~ 9^' °22 " *' 

^l2 l2 

we select the one that is consistent with the solution obtained from: 



that is, 

4i 



< 


= mj -^ t'2 


t\- 


:^0 




d'o 


= 1 




4i 


= - 


.'"112 
t'2 





^'13 = <3<^tr 






m'u 
m,, 


m,, 


^22 = 


'^22 

m'l, 
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(3) In most cases, the multiple ambiguity in the transformation angle, if it exists, can 
be resolved by the extra constraint equation; 

mj2 + ^21 = 0. 

In general, this constraint will be satisfied only for the correct rotation. We restate 
that these ambiguities occur only when the gaussian curvature at O is non-positive 
(the fixation point is either a parabolic cylinder or a hyperbolic paraboloid point). 
For surfaces with positive gaussian curvature at O, the solution for the transfor- 
mation angle is always unique. If we have not resolved the angle ambiguity yet, 
we, once again, turn to: 

since these constraints should only be satisfied for the correct transformation. 
It is very unlikely that the incorrect transformation satisfisfies all of the above 
constraints. 

(4) We then compute the motion and surface parameters in the original coordinate 
frame using the transformation equations given by: 

Wi = Qriw[, ti = flj/t', di = ?lud\, dim = ^il^jmdij- 



7. Summary 

In this paper, we have investigated the problem of recovering the motion of an observer 
relative to a stationary environment directly from the time-varying imagery. We assume 
that, in a small neighborhood, the surface of the scene can be approximated by a quadratic 
patch. We first compute 15 intermediate parameters termed essential parameters from 
a set of linear constraint equations in terms of the image intensity derivatives. We then 
solve for the observer motion parameters and the local structure of the scene (slope and 
curvature parameters) in closed form. 

We have shown that the solution is unique in every case. The previous results of 
Longuet-Higgins and Prazdny [10], or that of Waxman et al. [23,25] fall short of showing 
uniqueness for every case. Furthermore, they use the optical flow data. So far, no method 
for computing a correct optical flow is proposed in the literature than can be reliably used 
for the optical flow based methods. 

In a paper in press [17], we consider smooth textured surfaces with continuous n^ 
order derivatives (no restriction on n) . We show that the problem of recovering the motion 
and the coefficients of the surface can be decoupled into two problems: (1) determining 
the motion and the first and second order parameters of the surface through the procedure 
presented in this paper, and (2) using the solution from the first problem to solve for the 
higher order coefficients from a set of linear constraints. This suggests that including the 
higher order terms of the surface is not necessary for recovering the motion parameters. 
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Appendix 1: 

Earlier, we presented a procedure of recovering the motion of an observer relative to a 
locally quadratic patch. It is based on solving a set of nonlinear equations in a transformed 
image coordinate frame. To compute the transformation angle, we need to solve the 
following cubic equation: 

(^2tan^ - ii)(<fiitan^0 + 2di2tan^ + <f22) = 0. 

Only one of the solutions of this equation, tan^ = <i/<2, corresponds to the proper 
rotation. The other two solutions: 

-di2 ± Y<^12 ~ <^llC'22 

tan B — , 

an 

if they exist, are directly related to the local properties of the surface, namely the gaussian 

curvature of the surface, at (fixation point). Note that we are only interested in the 

real solutions of the above quadratic solution, and these may not exist. Here, we will 

investigate these in more details. 

As before, let Vi = (tand 1), and let: 

dXl dXidX2 



Hij — -dtj — 



d^Z d^Z 

\dX1dX2 dXl J 



represent the Hessian of the surface at O. The sign of the guassian curvature of the 
surface at this point is that of \Hij\, where: 

\Hij\ = dnd22 - di2- 

If \Hij\ > 0, that is, the surface has a positive guassian curvature at O (see Figure 
2a), then it is clear that the two solutions of 6: 

-di2 ± \ld\2 - dnd22 
tan 6 = ^4 

Oil 

are complex. The fixation point O, in this case, is referred to as a eliptic paraboloid point 
in the literature. 

Now consider the case when Hij is singular, that is \Hij\ = 0, but dij 7^ 0. In this 
case, the guassian curvature, at O is zero (see Figure 2b). In this case (O is a parabolic 
cylinder point), the two real solutions of 9 become identical, given by: 

^ <^12 <^22 

tan = — -r— = - -7— . 
dii di2 
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Gincave Surftoe 




Convex Surface 



Eliptic 

Paraboloid 
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(a) Pondve Gaussian Curvature at P 



Parabolic Qlinder 




(b) Zero Gaussain Curvature at P 




(c) Negative Gaussian Curvature at P 



Figure 2: Examples of Surfaces with Positive. Zero, and Negative Gaussian Curvature 



Finally consider the case when the surface has a negative gaussian curvature, that is 
\Hij\ < at O (see Figure 2c). The point of fixation O is referred to as a hyperbolic 
paraboloid or a saddle point. In this case, both solutions of 6 are real, and given by: 



tan^ 



■di2 ± \jd\2 - ^11^22 



d\\ 
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Appendix 2: 

In this appendix, we present the results of section 2 in the vector form for those readers 
that are not familiar with the tensor notation. 

2.1. Perspective Projection 

Let X denote the perspective projection of the point X onto the image plane (plane 2 = 1). 
Then: 

X • z 

2.2. Rigid Body Motion 

Let t and w denote the translational and rotational velocities of the observer relative to 

the scene, then: 

dX 

-r- = -w X X-t. 
at 

2.3. Image Motion Field and Optical Flow 

The projection of the 3-D velocity field, induced by a rigid body motion, onto the image 
is known as the image motion field, and is given by: 

dx d 1 
dt dfK • z 

For convenience, we introduce the notation Xt and X^ for the time derivatives of x and 
X, respectively. We then have 

'''^3ri^'"(xT^(^''^^^' 
which can also be written in the compact form: 

since a X (b X c) = (c • a)b - (a • b)c. Note that X = (X • z)x, and therefore, 

x« = :^jr-T(z X {Xt x x)). 
Finally, substituting for Xt, we obtain: 



Xt = -(z X [x X (x X (jj - :^—7t)) j . 



X z 

2.4. Brightness Change Constraint Equation 

The brightness of the image of a particular patch of a surface depends on many factors. 
It may for example vary with the orientation of the patch relative to the viewer or the 
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light source. In many cases, however, it remains at least approximately constant as the 
surface moves in the environment. If we assume that the image brightness of a patch 

remains constant, we have: 

dE 

— =0, 
dt 

or 

dE dx dE 

dx dt dt 

is the image brightness gradient. It is convenient to use the notation E'x for this quantity 

and Et for the time derivative of the brightness. Then we can write the brightness change 

equation in the simple form 

E^-yit + Et = 0. 

Substituting for X( we get 

Et- E^- (zx {xx {xxu- ^ry^))) = 0- 

If can be eaasily verified that if we define s = [E^ x z) x x, and v = -s x x, then the 
above equation can be written in the simple form: 

Et + y-u) + r^-rs -1=0. 
X • z 

This is the brightness change constraint equation in the case of rigid body motion. 



Appendix 3: 

We defined d=l/Z. Therefore: 

dd _ J_dZ^ d^d _ 1 d^Z 2 dZ dZ 

d^i~~Z^dxi' dxidxj ~ Z2 dxidxj '^ Z^^dx,'^ dxj ' ' 

Evaluating these expressions at the origin: 

^dxr' ~ Z2^dx/"' ^dx.dx/" " Z^^dx.dx/" ^ Z^^dx,'"^dxj'" 
We now use the perspective projection equation Xj — XiZ: 

dZ^ _ ^Z_dXrn _ ^Z_ 2 X —\ 

dxi dXm dxi dXm '"^ ^ dxi 

Therefore: 

dZ _ dZ __ dZ 1 

dx^ dXi dXm 
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Now: 

d^Z d , dZ ,, dZ. dZ ,, dZ ^ dZ d dZ . 

or, 

d'Z d ,dZ.dXn,, ^ dZ. dZ .^ dZ dZ d-'Z 

This simplifies to: 
d^Z d^Z , dZ,,^ ^ dZ , dZ ., dZ - dZ a^z 

After substituting for dZ/dx^, we arrive at: 
d^Z d^Z - - dZ r , dZ ,_.,, dZ , dZ . d'^Z - 

dXidXj dXmdXn ^ dXm O^k "Xj OXt OXiXj 

where, 

dZ , dZ ,_i 

Ztm^ kmZ + XmZ-^[l- Xk-^) ■ 

Solving for d'^Z/dxiXj, we obtain: 
d^Z dZ,_i- ~ d'^Z „. ^^x-2/c ^^/: ^\1^ 

ai^a;- ^ (' - ^^ax;^ ^^-^^'^ax;;:ax: + ^^'"'^^ax;^^ ^'^-ax/^^^^ax.^ax^- 

Evaluating {dZ/dxi)„ and {d"^ Z/dxidxj)„ at the origin results in: 

.az _ dz d'z , d'z ^2Zf^)f— ) 

^a^^" ~ "^"ax^; ^ax,x/" ~ ^-^ax.axy^" '^^^"^ax>'^ax/'■• 
Finally, substituting for {dZ/dxi),, and {d"^ Z/dx^dx^),, into the equations for {dd/dx,)„ 
and {d^d/dxidxj)„ we arrive at: 

*ai/"~ z.}dx,'"' '-dx.ax,'" '■dx.dXj'" 
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